Clear[initialcrossing,initialmatrix,braidword,Quad, A, F, HigherOrder,kerbeforesub,keraftersub,exp2,poly,poly2,tr,knot,n,p1,p2]

knot={3,1};
f[{i_,n_}]:=If[n>0,ConstantArray[{i,i+1},{n}],ConstantArray[{i+1,i},{-n}]]
list=Flatten[f/@KnotData[knot,"BraidWord"]];
initialcrossing=firstcrossing[Take[list,{1,2}]];
initialmatrix=firstmatrix[Take[list,{1,2}]];
braidword=Sequence[Partition[Take[list,2-Length[list]],2]];
n=KnotData[knot,"BraidIndex"];
firstcrossing[{x_,y_}]:=If[x<y,(Subscript[HO, x,y]),(Subscript[INV, y,x])] ;
firstmatrix[{x_,y_}]:= If[x<y,(Subscript[Matrix, x,y]),(Subscript[Imatrix, y,x])];
alex=KnotData[knot,"AlexanderPolynomial"][t^2];
amp=KnotData[knot,"Amphichiral"];
rep2[i_,j_]:=Which[i==j==0,1,True,((i!j!Coefficient[exp2,(Subscript[d, 2])^i (Subscript[z, 2])^j ])/.{(Subscript[z, 2])->0,(Subscript[d, 2])->0})]
rep3[i_,j_,k_,l_]:=Which[i==j==k==l==0,1,True,((i!j!k!l!Coefficient[exp2,(Subscript[d, 2])^i (Subscript[z, 2])^j (Subscript[d, 3 ])^k (Subscript[z, 3])^l])/.{(Subscript[z, 2])->0,(Subscript[d, 2])->0,(Subscript[z, 3])->0,(Subscript[d, 3])->0})]
rep4[i_,j_,k_,l_,m_,n_]:=Which[i==j==k==l==m==n==0,1,True,((i!j!k!l!m!n!Coefficient[exp2,(Subscript[d, 2])^i (Subscript[z, 2])^j (Subscript[d, 3 ])^k (Subscript[z, 3])^l (Subscript[d, 4 ])^m (Subscript[z, 4])^n])/.{(Subscript[z, 2])->0,(Subscript[d, 2])->0,(Subscript[z, 3])->0,(Subscript[d, 3])->0,(Subscript[z, 4])->0,(Subscript[d, 4])->0})]
rep5[i_,j_,k_,l_,m_,n_,o_,p_]:=Which[i==j==k==l==m==n==o==p==0,1,True,((i!j!k!l!m!n!o!p!Coefficient[exp2,(Subscript[d, 2])^i (Subscript[z, 2])^j (Subscript[d, 3 ])^k (Subscript[z, 3])^l (Subscript[d, 4 ])^m (Subscript[z, 4])^n (Subscript[d, 5 ])^o (Subscript[z, 5])^p])/.{(Subscript[z, 2])->0,(Subscript[d, 2])->0,(Subscript[z, 3])->0,(Subscript[d, 3])->0,(Subscript[z, 4])->0,(Subscript[d, 4])->0,(Subscript[z, 5])->0,(Subscript[d, 5])->0})]
rep6[i_,j_,k_,l_,m_,n_,o_,p_,q_,r_]:=Which[i==j==k==l==m==n==o==p==q==r==0,1,True,((i!j!k!l!m!n!o!p!q!r!Coefficient[exp2,(Subscript[d, 2])^i (Subscript[z, 2])^j (Subscript[d, 3 ])^k (Subscript[z, 3])^l (Subscript[d, 4 ])^m (Subscript[z, 4])^n (Subscript[d, 5 ])^o (Subscript[z, 5])^p (Subscript[d, 6])^q (Subscript[z, 6])^r])/.{(Subscript[z, 2])->0,(Subscript[d, 2])->0,(Subscript[z, 3])->0,(Subscript[d, 3])->0,(Subscript[z, 4])->0,(Subscript[d, 4])->0,(Subscript[z, 5])->0,(Subscript[d, 5])->0,(Subscript[z, 6])->0,(Subscript[d, 6])->0})]
(Subscript[trace, 2])[s_]:=Expand[Plus@@(s/.Rule[{a_,b_},c_]:>rep2[a,b]c)]
(Subscript[trace, 3])[s_]:=Expand[Plus@@(s/.Rule[{a_,b_,e_,f_},c_]:>rep3[a,b,e,f]c)]
(Subscript[trace, 4])[s_]:=Expand[Plus@@(s/.Rule[{a_,b_,e_,f_,g_,h_},c_]:>rep4[a,b,e,f,g,h]c)]
(Subscript[trace, 5])[s_]:=Expand[Plus@@(s/.Rule[{a_,b_,e_,f_,g_,h_,i_,j_},c_]:>rep5[a,b,e,f,g,h,i,j]c)]
(Subscript[trace, 6])[s_]:=Expand[Plus@@(s/.Rule[{a_,b_,e_,f_,g_,h_,i_,j_,k_,l_},c_]:>rep6[a,b,e,f,g,h,i,j,k,l]c)]
qh[n_]:=(1+h Coefficient[Product[1+2h (Subscript[z, i])(Subscript[d, i]),{i,2,n}],h,1]+h^2 Coefficient[Product[1+2h (Subscript[z, i])(Subscript[d, i])+2 h^2((Subscript[z, i])(Subscript[d, i])+(Subscript[z, i])^2(Subscript[d, i])^2),{i,2,n}],h,2]+h^3 Coefficient[Product[1+2h (Subscript[z, i])(Subscript[d, i])+2 h^2((Subscript[z, i])(Subscript[d, i])+(Subscript[z, i])^2(Subscript[d, i])^2)+(4/3)h^3((Subscript[z, i])(Subscript[d, i])+3((Subscript[z, i])^2)((Subscript[d, i])^2)+((Subscript[z, i])^3)(Subscript[d, i])^3),{i,2,n}],h,3]+h^4 Coefficient[Product[1+2h (Subscript[z, i])(Subscript[d, i])+2 h^2((Subscript[z, i])(Subscript[d, i])+(Subscript[z, i])^2(Subscript[d, i])^2)+(4/3)h^3((Subscript[z, i])(Subscript[d, i])+3((Subscript[z, i])^2)((Subscript[d, i])^2)+((Subscript[z, i])^3)(Subscript[d, i])^3)+(2/3)h^4((Subscript[z, i])(Subscript[d, i])+7((Subscript[z, i])^2)((Subscript[d, i])^2)+6((Subscript[z, i])^3)((Subscript[d, i])^3)+((Subscript[z, i])^4)((Subscript[d, i])^4)),{i,2,n}],h,4])
(Subscript[p, l_,k_])[i_,j_]:=Which[i==k&&j==l,t,i==l&&j==k,t,i==j==l,0,i==j==k,1-t^2,i==j,1,True,0];
(Subscript[Matrix, l_,k_]):=Array[(Subscript[p, l,k]),{n,n}];
(Subscript[Imatrix, l_,k_]):=Inverse[(Subscript[Matrix, l,k])];
(Subscript[Sz, l_,k_])[(Subscript[z, j_])]:=Expand[Sum[(Subscript[z, i])(Subscript[Imatrix, l,k])[[i,j]],{i,1,n}]];
(Subscript[Sder, l_,k_])[(Subscript[d, j_])]:=Expand[Sum[(Subscript[Matrix, l,k])[[i,j]](Subscript[d, i]),{i,1,n}]];
(Subscript[ISz, l_,k_])[(Subscript[z, j_])]:=Expand[Sum[(Subscript[z, i])(Subscript[Matrix, l,k])[[i,j]],{i,1,n}]];
(Subscript[ISder, l_,k_])[(Subscript[d, j_])]:=Expand[Sum[(Subscript[Imatrix, l,k])[[i,j]](Subscript[d, i]),{i,1,n}]];
Higher[a_,b_,w_,x_]:=1+h (2 a b w x-((2 (a^2) w x)/t)+2 (a^2) t w x+3 (a^2) (x^2)-(((a^2) (x^2))/(t^2))+((a b (x^2))/t)-3 a b t (x^2)-2 (a^2) (t^2) (x^2))+(h^2) (2 a b w x-((2 (a^2) w x)/t)+2 (a^2) t w x+2 (a^2) b (w^2) x-((2 (a^3) (w^2) x)/t)+2 (a^3) t (w^2) x+5 (a^2) (x^2)-(((a^2) (x^2))/(t^2))+((a b (x^2))/t)-5 a b t (x^2)-4 (a^2) (t^2) (x^2)+4 (a^3) w (x^2)+2 a (b^2) w (x^2)-((2 (a^2) b w (x^2))/t)-2 (a^2) b t w (x^2)-4 (a^3) (t^2) w (x^2)-4 (a^4) (w^2) (x^2)+2 (a^2) (b^2) (w^2) (x^2)+((2 (a^4) (w^2) (x^2))/(t^2))-((4 (a^3) b (w^2) (x^2))/t)+4 (a^3) b t (w^2) (x^2)+2 (a^4) (t^2) (w^2) (x^2)+(10/3) (a^2) b (x^3)-((2 (a^2) b (x^3))/(3 (t^2)))+((4 (a^3) (x^3))/(3 t))+((2 a (b^2) (x^3))/(3 t))-(16/3) (a^3) t (x^3)-(14/3) a (b^2) t (x^3)+(4/3) (a^2) b (t^2) (x^3)+4 (a^3) (t^3) (x^3)+14 (a^3) b w (x^3)+((2 (a^4) w (x^3))/(t^3))-((4 (a^3) b w (x^3))/(t^2))-((8 (a^4) w (x^3))/t)+((2 (a^2) (b^2) w (x^3))/t)+10 (a^4) t w (x^3)-6 (a^2) (b^2) t w (x^3)-10 (a^3) b (t^2) w (x^3)-4 (a^4) (t^3) w (x^3)+((13 (a^4) (x^4))/2)-3 (a^2) (b^2) (x^4)+(((a^4) (x^4))/(2 (t^4)))-(((a^3) b (x^4))/(t^3))-((3 (a^4) (x^4))/(t^2))+(((a^2) (b^2) (x^4))/(2 (t^2)))+((6 (a^3) b (x^4))/t)-11 (a^3) b t (x^4)-6 (a^4) (t^2) (x^4)+(9/2) (a^2) (b^2) (t^2) (x^4)+6 (a^3) b (t^3) (x^4)+2 (a^4) (t^4) (x^4))
Inv[a_,b_,w_,x_]:=1+h ((b^2) (w^2)+((a b (w^2))/t)+a b t (w^2)-(b^2) (t^2) (w^2)-2 a b w x)+(h^2) ((b^2) (w^2)-((a b (w^2))/t)+a b t (w^2)-(b^2) (t^2) (w^2)+(10/3) a (b^2) (w^3)+((4 a (b^2) (w^3))/(3 (t^2)))-((2 (a^2) b (w^3))/(3 t))+((8 (b^3) (w^3))/(3 t))+(2/3) (a^2) b t (w^3)-(8/3) (b^3) t (w^3)-(2/3) a (b^2) (t^2) (w^3)+(a^2) (b^2) (w^4)+(((b^4) (w^4))/2)+(((a^2) (b^2) (w^4))/(2 (t^2)))+((a (b^3) (w^4))/t)+(1/2) (a^2) (b^2) (t^2) (w^4)-(b^4) (t^2) (w^4)-a (b^3) (t^3) (w^4)+(1/2) (b^4) (t^4) (w^4)+2 a b w x+2 (a^2) b (w^2) x-4 (b^3) (w^2) x-((4 a (b^2) (w^2) x)/t)-4 a (b^2) t (w^2) x+4 (b^3) (t^2) (w^2) x-2 a (b^3) (w^3) x-((2 (a^2) (b^2) (w^3) x)/t)-2 (a^2) (b^2) t (w^3) x+2 a (b^3) (t^2) (w^3) x+2 a (b^2) w (x^2)+2 (a^2) (b^2) (w^2) (x^2))
(Subscript[HO, i_,j_]):=Higher[(Subscript[z, i]),(Subscript[z, j]),(Subscript[d, i]),(Subscript[d, j])]
(Subscript[INV, i_,j_]):=Inv[(Subscript[z, i]),(Subscript[z, j]),(Subscript[d, i]),(Subscript[d, j])]
(Subscript[Der, i_])[f_]:=f+Sum[(1/k!)D[f,{(Subscript[w, i]),k},{(Subscript[x, i]),k}],{k,1,8}]
NO3[l_,{i_,j_}]/;i<j :=((Subscript[Der, j])[(Subscript[Der, i])[Expand[1+h Coefficient[((l/.Flatten[Table[{(Subscript[z, k])->(Subscript[Sz, i,j])[(Subscript[z, k])],(Subscript[d, k])->(Subscript[Sder, i,j])[(Subscript[d, k])]},{k,1,n}]])/.{(Subscript[d, i])->(Subscript[w, i]),(Subscript[d, j])->(Subscript[w, j])})*Higher[(Subscript[x, i]),(Subscript[x, j]),(Subscript[d, i]),(Subscript[d, j])],h,1]+h^2Coefficient[((l/.Flatten[Table[{(Subscript[z, k])->(Subscript[Sz, i,j])[(Subscript[z, k])],(Subscript[d, k])->(Subscript[Sder, i,j])[(Subscript[d, k])]},{k,1,n}]])/.{(Subscript[d, i])->(Subscript[w, i]),(Subscript[d, j])->(Subscript[w, j])})*Higher[(Subscript[x, i]),(Subscript[x, j]),(Subscript[d, i]),(Subscript[d, j])],h,2]]]])/.{(Subscript[w, i])->(Subscript[d, i]),(Subscript[w, j])->(Subscript[d, j]),(Subscript[x, i])->(Subscript[z, i]),(Subscript[x, j])->(Subscript[z, j])}
NO3[l_,{i_,j_}]/;i>j:=((Subscript[Der, j])[(Subscript[Der, i])[Expand[1+h Coefficient[((l/.Flatten[Table[{(Subscript[z, k])->(Subscript[ISz, j,i])[(Subscript[z, k])],(Subscript[d, k])->(Subscript[ISder, j,i])[(Subscript[d, k])]},{k,1,n}]])/.{(Subscript[d, i])->(Subscript[w, i]),(Subscript[d, j])->(Subscript[w, j])})*Inv[(Subscript[x, j]),(Subscript[x, i]),(Subscript[d, j]),(Subscript[d, i])],h,1]+h^2Coefficient[((l/.Flatten[Table[{(Subscript[z, k])->(Subscript[ISz, j,i])[(Subscript[z, k])],(Subscript[d, k])->(Subscript[ISder, j,i])[(Subscript[d, k])]},{k,1,n}]])/.{(Subscript[d, i])->(Subscript[w, i]),(Subscript[d, j])->(Subscript[w, j])})*Inv[(Subscript[x, j]),(Subscript[x, i]),(Subscript[d, j]),(Subscript[d, i])],h,2]]]])/.{(Subscript[w, i])->(Subscript[d, i]),(Subscript[w, j])->(Subscript[d, j]),(Subscript[x, i])->(Subscript[z, i]),(Subscript[x, j])->(Subscript[z, j])}

subd=Flatten[Table[(Subscript[d, i])->(Subscript[w, i]),{i,2,n}]];
subz=Flatten[Table[(Subscript[z, i])->(Subscript[x, i]),{i,2,n}]];

g[l_,{x_ ,y_}]:=If[x<y,l.((Subscript[Matrix, x,y])),l.((Subscript[Imatrix, y,x]))]

qhend3[f_]:=Expand[(f/.subd)+h (Coefficient[qh[n],h,1]/.subz)+h^2 (Coefficient[qh[n],h,2]/.subz)+h^2 ((Coefficient[f,h,1]/.subd)*(Coefficient[qh[n],h,1]/.subz))]

Quad=Expand[Fold[g,initialmatrix,braidword]];
(Subscript[a, i_,j_]):=If[i==j,Quad[[i,j]]-1,Quad[[i,j]]];
A=Table[(Subscript[a, i,j]),{i,2,n},{j,2,n}];
F=Inverse[-A];

AbsoluteTiming[HigherOrder3=Fold[NO3,initialcrossing,braidword];]

AbsoluteTiming[kerbeforesub3=((l=qhend3[HigherOrder3];Do[l=(Subscript[Der, i])[l],{i,2,n,1}];l)/.Flatten[Table[{(Subscript[x, i])->(Subscript[z, i]),(Subscript[w, i])->(Subscript[d, i])},{i,2,n}]]);]

sub=Flatten[Table[(Subscript[z, k])->Expand[Sum[Quad[[i,k]] (Subscript[z, i]),{i,1,n}]],{k,1,n}]]

AbsoluteTiming[keraftersub=1+h Expand[Coefficient[kerbeforesub3/.sub,h,1]]+h^2 Expand[Coefficient[kerbeforesub3/.sub,h,2]];]

sub3=Flatten[Table[{(Subscript[z, i])->(Subscript[z, i])+Expand[Sum[F[[k,i-1]]Quad[[1,k+1]](Subscript[z, 1]),{k,1,n-1}]],(Subscript[d, i])->(Subscript[d, i])+Expand[Sum[F[[i-1,l]]Quad[[l+1,1]](Subscript[d, 1]),{l,1,n-1}]]},{i,2,n}]]

AbsoluteTiming[twopoly=1+h Expand[Coefficient[keraftersub/.sub3,h,1]]+h^2 Expand[Coefficient[keraftersub/.sub3,h,2]];]

sub4=Sum[(Expand[Sum[F[[i-1,j-1]](Subscript[d, i])(Subscript[z , j]),{j,2,n},{i,2,n}]])^i/i!,{i,1,10}]

AbsoluteTiming[exp2=Expand[1+sub4];]

sub5=Flatten[Table[{(Subscript[d, i]),(Subscript[z, i])},{i,2,n}]]

AbsoluteTiming[poly2=CoefficientRules[twopoly,sub5]]

AbsoluteTiming[tr=(Subscript[trace, n])[poly2];]

test1=Together[Coefficient[tr,h,1]]

test2=Together[Coefficient[tr,h,2]]

p=If[amp==True,Cancel[Denominator[test1]/alex],Cancel[Denominator[test1]/alex^2]]
q=If[amp==True,Cancel[Denominator[test2]/alex^3],Cancel[Denominator[test2]/alex^4]]

series=If[amp==True,Normal[Series[((((1+h test1+h^2test2+h^3 test3+h^4 test4)(1/alex))/.{t->t Exp[h]})(alex)/.{h-> (h/2)-((h^2)/4)})/.{h->-Sqrt[1+h]+1/Sqrt[1+h]},{h,0,2}]],Normal[Series[((((1+h test1+h^2test2+h^3 test3+h^4 test4)(1/alex))/.{t->t Exp[h]})alex)/.{h-> (h/2)-((h^2)/4)},{h,0,2}]]]

p1=Expand[Numerator[Coefficient[series,h,1]]/p]

p2=Expand[Numerator[Coefficient[series,h,2]]/q]

SessionTime[]